Check residuals

cohortfile <- read.csv("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/cohortfiles/dti_fa/dti_fa_cohortfile.csv")
cohortfile <- cohortfile %>% rename(subjectID=sub) %>% select(subjectID, age, sex, mean_fd)

tract_profiles <- fread("/cbica/projects/luo_wm_dev/input/HCPD/HCPD_tractprofiles/all_subjects/collated_tract_profiles_reoriented.tsv")

cohortfile <- cohortfile %>% select(subjectID, age, sex, mean_fd)
gam_df <- merge(tract_profiles, cohortfile, by="subjectID")
gam_df <- gam_df %>% mutate(tract_node = paste0(tract_hemi, "_", nodeID))
gam_df$sex <- as.factor(gam_df$sex)

subjects <- read.table("/cbica/projects/luo_wm_dev/input/HCPD/subject_list/HCPD_subject_list_N569.txt")
subjects <- subjects %>% setNames("subjectID")
subject_subset <- subjects$subjectID[c(1:75)]
gam_df <- gam_df %>% filter(subjectID %in% subject_subset) # fit GAMs on a subset of 75 participants

fit_tractwise_gamm_re <- function(tract_name, scalar) {
  df <- gam_df %>% filter(tract_hemi == tract_name)
  df$subjectID <- as.factor(df$subjectID)
  te_gam <- gamm(as.formula(sprintf("%s ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') + 
                                  ti(age, nodeID, k=c(3, 24), fx = F) + sex + mean_fd", scalar)), random = list(subjectID=~1), data = df) 
  
  print(summary(te_gam$gam))
  print(k.check(te_gam$gam))
  print(paste(tract_name, "done"))
 
  
  return(te_gam)  # return the model object 
}


fit_tractwise_gam_re <- function(tract_name, scalar) {
  df <- gam_df %>% filter(tract_hemi == tract_name)
  df$subjectID <- as.factor(df$subjectID)
  te_gam <- gam(as.formula(sprintf("%s ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') + 
                                  ti(age, nodeID, k=c(3, 24), fx = F) + sex + mean_fd + s(subjectID, bs = 're')", scalar)), data = df) 
  
  print(summary(te_gam))
  print(k.check(te_gam))
  print(paste(tract_name, "done"))
  
  
  return(te_gam)  # return the model object 
}

fit_tractwise_bam <- function(tract_name, scalar) {
  df <- gam_df %>% filter(tract_hemi == tract_name)
  df$subjectID <- as.factor(df$subjectID)
  te_gam <- bam(as.formula(sprintf("%s ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') + 
                                  ti(age, nodeID, k=c(3, 24), fx = F) +  s(nodeID, by = subjectID) + sex + mean_fd", scalar)), data = df, method = "fREML", discrete = TRUE) 
  
  print(summary(te_gam))
  print(k.check(te_gam))
  print(paste(tract_name, "done"))
 
  
  return(te_gam)  # return the model object 
}

# 3) Using bam and using separate smooth function for subjectID
#`bam(dti_md ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') + ti(age, nodeID, k=c(3, 24), fx = F) + sex + mean_fd + s(nodeID, by = subjectID), data = df, method = "fREML", discrete = TRUE)`

fit_tractwise_bam2 <- function(tract_name, scalar) {
  df <- gam_df %>% filter(tract_hemi == tract_name)
  df$subjectID <- as.factor(df$subjectID)
  te_gam <- bam(as.formula(sprintf("%s ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr', by = subjectID) + 
                                  ti(age, nodeID, k=c(3, 24), fx = F) + sex + mean_fd", scalar)), data = df, method = "fREML", discrete = TRUE) 
  
  print(summary(te_gam))
  print(k.check(te_gam))
  print(paste(tract_name, "done"))
 
  
  return(te_gam)  # return the model object 
}

tractwise_gamm_re <- lapply(unique(gam_df$tract_hemi), fit_tractwise_gamm_re, scalar="dti_md")
names(tractwise_gamm_re) <- unique(gam_df$tract_hemi)
saveRDS(tractwise_gamm_re, "/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/check_resid_tractwise_gamm_re.RData")

tractwise_gam_re <- lapply(unique(gam_df$tract_hemi), fit_tractwise_gam_re, scalar="dti_md")
names(tractwise_gam_re) <- unique(gam_df$tract_hemi)
saveRDS(tractwise_gam_re, "/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/check_resid_tractwise_gam_re.RData")

tractwise_bam <- lapply(unique(gam_df$tract_hemi), fit_tractwise_bam, scalar="dti_md")
names(tractwise_bam) <- unique(gam_df$tract_hemi)
saveRDS(tractwise_bam, "/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/check_resid_tractwise_bam.RData")


tractwise_bam2 <- lapply(unique(gam_df$tract_hemi), fit_tractwise_bam2, scalar="dti_md")
names(tractwise_bam2) <- unique(gam_df$tract_hemi)
saveRDS(tractwise_bam2, "/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/check_resid_tractwise_bam2.RData")

tractwise_gamm_re <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/check_resid_tractwise_gamm_re.RData")
tractwise_gam_re <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/check_resid_tractwise_gam_re.RData")
#tractwise_bam <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/check_resid_tractwise_gam_smooth.RData")

tractwise_bam <- readRDS("/cbica/projects/luo_wm_dev/output/HCPD/tract_profiles/GAM/dti_md/check_resid_tractwise_bam2.RData")

I fit 3 tractwise models on a subset of HCP-D (full sample = 569; subsetted sample = 75 that represents the full age range)

1. Using gam with subjectID as a random effect: gam(dti_md ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') + ti(age, nodeID, k=c(3, 24), fx = F) + sex + mean_fd + s(subjectID, bs = 're'), data = df)

2. Using gamm and having subjectID as a random effect (gamm relies on lme, I believe): gamm(dti_md ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr') + ti(age, nodeID, k=c(3, 24), fx = F) + sex + mean_fd, random = list(subjectID=~1), data = df)

3. Using bam and using separate smooth function for subjectID but with k=24 bam(dti_md ~ s(age, k = 3, fx = T, bs = 'cr') + s(nodeID, k = 24, fx = F, bs = 'cr', by = subjectID) + ti(age, nodeID, k=c(3, 24), fx = F) + sex + mean_fd, data = df, method = "fREML", discrete = TRUE)

plot_compare_residuals <- function(tract, num_subjects) { 
  gam_df_toplot <- gam_df %>% filter(tract_hemi == tract)
  
  # Extract residuals for each model
  gam_df_toplot$residuals_gam <- resid(tractwise_gam_re[[tract]])
  gam_df_toplot$residuals_gamm <- resid(tractwise_gamm_re[[tract]]$gam)
  gam_df_toplot$residuals_bam <- resid(tractwise_bam[[tract]])
  #gam_df_toplot$residuals_bam2 <- resid(tractwise_bam2[[tract]])
  
  subjects_to_plot <- subject_subset[c(1:num_subjects)]
    
  # Plot residuals for each model
  gam_plot <- ggplot(gam_df_toplot[gam_df_toplot$subjectID %in% subjects_to_plot, ], aes(x = nodeID, y = residuals_gam, color = factor(subjectID))) +
    geom_line() + theme(plot.title = element_text(hjust=0.5, size = 16, face="bold"), 
                        axis.title.x = element_blank(), 
                        axis.title.y = element_blank()) + 
    labs(title = "Raw Residual Curves from Model 1 \ngam: s(subjectID, bs = 're')", x = "Node ID", y = "Residuals")  
  
  gamm_plot <- ggplot(gam_df_toplot[gam_df_toplot$subjectID %in% subjects_to_plot, ], aes(x = nodeID, y = residuals_gamm, color = factor(subjectID))) +
    geom_line() + theme(plot.title = element_text(hjust=0.5, size = 16, face="bold"),
                        axis.title.x = element_blank(), 
                        axis.title.y = element_blank()) + 
    labs(title = "Raw Residual Curves from Model 2 \ngamm: random = list(subjectID=~1)", x = "Node ID", y = "Residuals")  
  
  # Plot smooth residuals from the new approach
  bam_plot <- ggplot(gam_df_toplot[gam_df_toplot$subjectID %in% subjects_to_plot, ], aes(x = nodeID, y = residuals_bam, color = factor(subjectID))) +
    geom_line() + theme(plot.title = element_text(hjust=0.5, size = 16, face = "bold"),
                        axis.title.x = element_blank(), 
                        axis.title.y = element_blank()) + 
    labs(title = "Smooth Residual Fit from Model 3 \nbam: s(nodeID, k=24, by = subjectID)", x = "Node ID", y = "Residuals")  
  
  #bam2_plot <- ggplot(gam_df_toplot[gam_df_toplot$subjectID %in% subjects_to_plot, ], aes(x = nodeID, y = residuals_bam2, color = factor(subjectID))) +
   # geom_line() + theme(plot.title = element_text(hjust=0.5, size = 16, face = "bold"),
      #                  axis.title.x = element_blank(), 
    #                    axis.title.y = element_blank()) + 
    #labs(title = "Smooth Residual Fit from Model 4 \nbam: s(nodeID, k=24, by = subjectID)", x = "Node ID", y = "Residuals")
  
  x.grob <- textGrob("Position on Tract (Node ID)", 
                   gp=gpar(col="black", fontsize=16))
  y.grob <- textGrob("Residuals", 
                   gp=gpar(col="black", fontsize=16), rot=90)
  
  z.grob <- textGrob(tract, 
                   gp=gpar(col="black", fontsize=20))
  
  plots <- ggarrange(gam_plot, gamm_plot, bam_plot, ncol=3, common.legend=T, legend = "right") 
  


  return(grid.arrange(arrangeGrob(plots, left = y.grob, bottom = x.grob, top = z.grob)))
}

Residual plots for each type of model for each tract (showing 5 subjects only)

tracts <- unique(gam_df$tract_hemi)
tracts <- tracts[-c(which(tracts=="Forceps_Major"))]
lapply(tracts, plot_compare_residuals, 5)

## [[1]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[2]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[3]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[4]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[5]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[6]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[7]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[8]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[9]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[10]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[11]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[12]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[13]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[14]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[15]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[16]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[17]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[18]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[19]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[20]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[21]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]

Residual plots for each type of model for each tract (showing 15 subjects – idk this more helpful to look at?)

lapply(tracts, plot_compare_residuals, 15)

## [[1]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[2]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[3]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[4]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[5]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[6]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[7]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[8]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[9]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[10]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[11]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[12]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[13]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[14]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[15]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[16]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[17]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[18]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[19]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[20]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]
## 
## [[21]]
## TableGrob (1 x 1) "arrange": 1 grobs
##   z     cells    name            grob
## 1 1 (1-1,1-1) arrange gtable[arrange]

ACF plots for each type of model

plot_compare_acf <- function(tract) {
  # Define the layout: 2 rows (1 for title, 1 for plots), 3 columns for plots
layout(matrix(c(1, 1, 1, 2, 3, 4), ncol = 3, byrow = TRUE), heights = c(1, 4))

# Plot the title
par(mar = c(0, 0, 1, 0))  # Set margins for the title
plot.new()
text(0.5, 0.5, tract, cex = 2, font = 2)

# Reset margins and set up for 1 row, 3 columns for the plots
par(mar = c(5, 4, 4, 2) + 0.1)

# Plot the ACFs
acf(resid(tractwise_gam_re[[tract]], type = "response"), main = paste0("Model 1 \ngam: s(subjectID, bs = 're')"))
acf(resid(tractwise_gamm_re[[tract]]$gam, type = "response"), main = paste0("Model 2 \ngamm: random = list(subjectID=~1)"))
acf(resid(tractwise_bam[[tract]], type = "response"), main = paste0("Model 3 \nbam: s(nodeID, k=24, by = subjectID)"))
}
 

#plot_compare_acf4 <- function(tract) {
  # Define the layout: 2 rows (1 for title, 1 for plots), 3 columns for plots
#  layout(matrix(c(1, 1, 1, 1, 2, 3, 4, 5), ncol = 4, byrow = TRUE), heights = c(1, 4))
  
  # Plot the title
#  par(mar = c(0, 0, 1, 0))  # Set margins for the title
#  plot.new()
#  text(0.5, 0.5, tract, cex = 2, font = 2)
  
  # Reset margins and set up for 1 row, 4 columns for the plots
#  par(mar = c(5, 4, 4, 2) + 0.1)
  
  # Plot the ACFs
#  acf(resid(tractwise_gam_re[[tract]], type = "response"), main = paste0("Model 1 \ngam: s(subjectID, bs = 're')"))
#  acf(resid(tractwise_gamm_re[[tract]]$gam, type = "response"), main = paste0("Model 2 \ngamm: random = list(subjectID=~1)"))
#  acf(resid(tractwise_bam[[tract]], type = "response"), main = paste0("Model 3 \nbam: s(nodeID, by = subjectID)"))
  #acf(resid(tractwise_bam2[[tract]], type = "response"), main = paste0("Model 4 \nbam: s(nodeID, k=24, by = subjectID)"))

#}

lapply(unique(gam_df$tract_hemi), plot_compare_acf)

## [[1]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.978 0.948 0.915 0.881 0.849 0.820 0.797 0.779 0.767 0.759 0.755 0.753 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.751 0.750 0.749 0.747 0.743 0.739 0.732 0.724 0.714 0.702 0.689 0.676 0.662 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.649 0.637 0.626 0.616 0.607 0.599 0.590 0.581 0.572 0.561 0.551 0.539 0.528 
## 
## [[2]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.980 0.955 0.933 0.912 0.892 0.875 0.860 0.848 0.837 0.827 0.818 0.809 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.800 0.791 0.781 0.770 0.760 0.748 0.736 0.722 0.707 0.692 0.676 0.661 0.647 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.634 0.623 0.613 0.604 0.595 0.587 0.577 0.568 0.558 0.548 0.538 0.529 0.519 
## 
## [[3]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.972 0.942 0.919 0.897 0.873 0.851 0.831 0.815 0.802 0.792 0.783 0.776 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.770 0.763 0.757 0.749 0.740 0.730 0.719 0.707 0.695 0.682 0.670 0.658 0.647 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.636 0.626 0.617 0.607 0.599 0.590 0.582 0.573 0.565 0.556 0.547 0.537 0.527 
## 
## [[4]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.973 0.944 0.921 0.899 0.875 0.852 0.832 0.815 0.801 0.790 0.780 0.772 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.764 0.757 0.748 0.740 0.730 0.720 0.709 0.698 0.686 0.674 0.661 0.649 0.638 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.626 0.615 0.604 0.594 0.584 0.574 0.565 0.555 0.546 0.536 0.527 0.517 0.507 
## 
## [[5]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.627 0.688 0.618 0.550 0.512 0.504 0.514 0.531 0.548 0.563 0.567 0.561 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.548 0.534 0.519 0.505 0.495 0.489 0.485 0.480 0.475 0.469 0.463 0.456 0.447 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.438 0.430 0.423 0.417 0.412 0.406 0.399 0.390 0.382 0.374 0.367 0.363 0.360 
## 
## [[6]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.606 0.664 0.589 0.512 0.473 0.467 0.483 0.506 0.528 0.543 0.548 0.544 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.531 0.515 0.496 0.480 0.469 0.463 0.461 0.463 0.465 0.464 0.461 0.455 0.446 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.437 0.428 0.421 0.415 0.409 0.401 0.391 0.379 0.369 0.361 0.356 0.353 0.350 
## 
## [[7]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.964 0.908 0.855 0.809 0.778 0.763 0.761 0.765 0.772 0.776 0.776 0.770 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.760 0.747 0.732 0.716 0.699 0.683 0.667 0.651 0.636 0.623 0.612 0.607 0.608 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.613 0.617 0.616 0.609 0.595 0.577 0.557 0.539 0.522 0.509 0.500 0.494 0.491 
## 
## [[8]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.961 0.899 0.842 0.795 0.761 0.744 0.739 0.742 0.749 0.757 0.762 0.764 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.760 0.748 0.729 0.707 0.685 0.665 0.649 0.637 0.627 0.619 0.612 0.607 0.604 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.604 0.603 0.599 0.591 0.579 0.564 0.547 0.530 0.513 0.497 0.482 0.468 0.457 
## 
## [[9]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.964 0.917 0.880 0.852 0.829 0.808 0.790 0.775 0.763 0.754 0.748 0.745 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.743 0.742 0.740 0.737 0.731 0.723 0.713 0.701 0.689 0.677 0.665 0.654 0.643 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.633 0.625 0.618 0.611 0.606 0.601 0.595 0.589 0.583 0.576 0.568 0.559 0.550 
## 
## [[10]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.965 0.916 0.872 0.835 0.803 0.776 0.753 0.735 0.721 0.712 0.706 0.704 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.704 0.705 0.706 0.705 0.702 0.696 0.687 0.676 0.662 0.647 0.632 0.617 0.603 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.592 0.583 0.575 0.569 0.563 0.557 0.551 0.545 0.538 0.532 0.524 0.515 0.505 
## 
## [[11]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.963 0.930 0.913 0.899 0.881 0.861 0.842 0.825 0.811 0.799 0.789 0.780 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.772 0.763 0.756 0.747 0.738 0.729 0.720 0.710 0.701 0.691 0.681 0.671 0.661 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.651 0.642 0.633 0.625 0.617 0.610 0.603 0.596 0.589 0.582 0.575 0.568 0.560 
## 
## [[12]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.963 0.931 0.912 0.894 0.875 0.855 0.837 0.824 0.814 0.806 0.800 0.795 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.789 0.784 0.777 0.769 0.760 0.750 0.739 0.727 0.715 0.703 0.691 0.680 0.669 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.659 0.649 0.640 0.631 0.622 0.614 0.606 0.597 0.589 0.580 0.571 0.561 0.551 
## 
## [[13]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.964 0.936 0.916 0.894 0.871 0.851 0.836 0.825 0.816 0.809 0.801 0.793 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.785 0.776 0.766 0.755 0.745 0.734 0.723 0.712 0.702 0.693 0.684 0.675 0.667 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.659 0.650 0.642 0.633 0.624 0.614 0.604 0.594 0.584 0.575 0.567 0.558 0.549 
## 
## [[14]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.957 0.921 0.893 0.865 0.838 0.818 0.806 0.799 0.795 0.791 0.787 0.782 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.777 0.769 0.760 0.750 0.739 0.728 0.717 0.707 0.696 0.686 0.677 0.667 0.658 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.648 0.639 0.629 0.619 0.609 0.598 0.587 0.576 0.565 0.555 0.545 0.536 0.526 
## 
## [[15]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.966 0.934 0.909 0.889 0.872 0.855 0.838 0.820 0.804 0.790 0.777 0.765 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.755 0.746 0.737 0.729 0.722 0.714 0.707 0.700 0.693 0.685 0.677 0.667 0.658 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.647 0.637 0.626 0.616 0.605 0.594 0.582 0.571 0.560 0.548 0.537 0.526 0.515 
## 
## [[16]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.968 0.937 0.912 0.892 0.875 0.859 0.844 0.831 0.818 0.807 0.798 0.789 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.781 0.773 0.766 0.760 0.753 0.747 0.741 0.735 0.728 0.721 0.714 0.707 0.699 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.691 0.682 0.674 0.665 0.655 0.646 0.636 0.626 0.615 0.605 0.594 0.585 0.575 
## 
## [[17]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.974 0.948 0.928 0.906 0.881 0.856 0.834 0.816 0.801 0.788 0.778 0.769 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.762 0.756 0.750 0.746 0.741 0.736 0.730 0.723 0.714 0.703 0.691 0.678 0.664 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.650 0.636 0.622 0.609 0.597 0.586 0.576 0.567 0.559 0.551 0.543 0.535 0.526 
## 
## [[18]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.938 0.834 0.741 0.682 0.661 0.666 0.682 0.697 0.703 0.702 0.696 0.692 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.693 0.698 0.702 0.697 0.683 0.661 0.636 0.614 0.598 0.587 0.580 0.575 0.570 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.564 0.558 0.549 0.541 0.538 0.538 0.539 0.536 0.529 0.515 0.499 0.481 0.465 
## 
## [[19]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.979 0.958 0.938 0.920 0.903 0.886 0.870 0.855 0.842 0.832 0.823 0.816 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.810 0.804 0.799 0.793 0.786 0.779 0.771 0.763 0.755 0.746 0.737 0.728 0.719 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.710 0.700 0.691 0.681 0.671 0.661 0.651 0.642 0.632 0.622 0.613 0.603 0.594 
## 
## [[20]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.980 0.959 0.939 0.921 0.903 0.887 0.871 0.857 0.845 0.834 0.825 0.816 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.809 0.801 0.793 0.785 0.776 0.767 0.758 0.749 0.739 0.730 0.720 0.711 0.701 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.691 0.682 0.672 0.663 0.654 0.645 0.638 0.630 0.622 0.615 0.607 0.599 0.592 
## 
## [[21]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.977 0.953 0.930 0.908 0.888 0.869 0.851 0.835 0.821 0.808 0.797 0.787 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.778 0.770 0.762 0.755 0.748 0.741 0.734 0.726 0.718 0.710 0.701 0.692 0.682 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.673 0.663 0.653 0.642 0.632 0.622 0.612 0.602 0.593 0.584 0.575 0.566 0.556 
## 
## [[22]]
## 
## Autocorrelations of series 'resid(tractwise_bam[[tract]], type = "response")', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11    12 
## 1.000 0.987 0.971 0.957 0.944 0.931 0.918 0.905 0.892 0.879 0.866 0.854 0.843 
##    13    14    15    16    17    18    19    20    21    22    23    24    25 
## 0.832 0.822 0.812 0.803 0.794 0.784 0.775 0.766 0.756 0.747 0.737 0.728 0.720 
##    26    27    28    29    30    31    32    33    34    35    36    37    38 
## 0.711 0.703 0.694 0.685 0.676 0.667 0.658 0.649 0.640 0.632 0.624 0.617 0.609